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ABSTRACT 

We present radiation transfer simulations of a massive (SM©) protostar forming from a massive 
(Mc = 6OM0) protostellar core, extending the model developed by Zhang & Tan (2011, Paper I). The 
two principal improvements are (1) developing a model for the density and velocity structure of a disk 
wind that fills the bipolar outflow cavities, based in part on the disk- wind model of Blandford fc: Pavnd 
(|1982l ): and (2) solving for the radially varying accretion rate in the disk due to a supply of mass and 
angular momentum from the infall envelope and their loss to the disk wind. One consequence of the 
launching of the disk wind is a reduction in the amount of accretion power that is radiated by the disk. 
We also include a non-Keplerian potential appropriate for a growing, massive disk. For the transition 
from dusty to dust-free conditions where gas opacities dominate, we now implement a gradual change 
as a more realistic approximation of dust destruction. We study how the above effects, especially the 
outflow, influence the spectral energy distributions (SEDs) and the synthetic images of the protostar. 
Dust in the outflow cavity signiflcantly affects the SEDs at most viewing angles. It further attenuates 
the short-wavelength flux from the protostar, controlling how the accretion disk may be viewed, and 
contributes a signiflcant part of the near- and mid-IR fluxes. These fluxes warm the disk, boosting 
the mid- and far-IR emission. We flnd that for near face-on views, i.e. looking down the outflow 
cavity (although not too close to the axis), the SED from the near-IR to about 60 jiia is very flat, 
which may be used to identify such systems. We show that the near- facing outflow cavity and its walls 
are still the most signiflcant features in images up to 70 /im, dominating the mid-IR emission and 
determining its morphology. The thermal emission from the dusty outflow itself dominates the flux 
at 20/im. The detailed distribution of the dust in the outflow affects the morphology, so the model 
can be constrained by considering detailed intensity proflles along and perpendicular to the outflow 
axis. For example, even though the outflow cavity is wide, at 10 to 20 /im, the dust in the disk wind 
can make the outflow appear narrower than in the near-IR bands. 
Subject headings: ISM: clouds, dust, extinction, jets and outflows — stars: formation 



1. INTRODUCTION 

There is still no consensus on how massive stars form. 
Possible mechanisms include collapse from massive cores 
supported by turbulence and magnetic fields, as in the 
Turbulent Core model (McKee TanI I2003L hereafter 
MT03), competiti ve accretion at the centers of form- 
ing sta r clusters rBonnell et al.l l2001f ) , and stellar col- 
lisions (|Bonnell e t al. 1998^^ Theoretically, if a mas- 
sive star forms similarly to its low-mass counterpart, 
i.e. from a core with a typical ratio of rotational to 
gravi tational energy of a few percent ([Goodman et al.l 
1 19931 ). a massive accretion disk with a size of sev- 
eral hundreds to ~ 1000 AU and bipolar outfiows 
are expected to appear around the massive protostar 
|Tan 2005). Such disks and outflows have been seen 
in hydrodynamic simulations with varying implementa- 
tions of magnetic fleld, turbulence and radiation (e.g., 
iKrumholz etaIll2007L iKuiper et ani2010L iSeifried et atl 



[20Tll . [20T2l . [Peters et al.l[20TTL IVaidya et al.l[20TTf ). How- 
ever, resolving this issue observationally is diflicult as 
massive stars tend to form in distant, crowded, and 
highly obscured environments. Hot cores with collimated 
molecular outflows and disk-like structures have been 
found around massive young stars, but so far these ob- 
served disk-like structures (or toroids) are much larger 
and more massive than expected for the rotationally- 
supported disk of an individual protostar, and there is lit- 
tle dire ct evidence fo r Keplerian or near-Keplerian disks 
(e.g., Cesaroni et al.|[2Q0 7). The outflows have also been 
observed in near- and mid-IR continuum emission, and 
it has been argued that the out flow cavity may dominate 
the mid-IR morphology (.De Buizerll2006f ). 

In order to better understand the properties of these 
observed massive young stellar objects (MYSOs), a num- 
ber of radiative transfer models have been developed to 
compare with observations. One of the most popular 
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is the model grid by iRobitaihe et all (|2QQ6f ), but it is fo- 
cussed on low-mass star formation. Massive star-forming 
cores usually reside in clumps within giant molecular 
clouds with large mass surface densities, Sci ~ 1 g cm~^ 
(MT03), so they are themselves relatively dense and col- 
lapse with high accretion rates. The expected proper- 
ties of the fiducial c ores of MTQ3 l i e out side the pa- 
rameter space of the iRobitaille et al.l ()2QQ 6) grid. Also 
the components in this radiative transfer model are rela- 
tively simple — e.g., the density in the outflow is usually 
assumed to be either a constant value or a power-law. 
[Offner et^al] (j2Q12f ) compared this model grid with re- 
sults of numerical simulations of low-mass star forma- 
tion, and found possible discrepancies between the pa- 
rameters inferred from fitting with the radiative-transfer 
model grid and those of the protostars in the simulation, 
especially for deeply embedded objects. Ot her radiative- 
transf er model grids, such as those of Molinari et al.l 
()2QQ8f ). are for MYSOs, but the components are also rel- 
atively simple. So it is worth including a more realistic 
density distribution in the outflow cavity and construct- 
ing a fully self-consistent model for all the components 
of the star-forming core. 

This is the second of a series of papers on model- 
ing the spectral energy distributions (SEDs) and images 
of massive sta r-forming cores. In the previous paper 
(jZhang fc Tanl[2Ql l. hereafter Paper I), we studied a core 
with an initial mass of 60 Mq that contained an 8 Mq 
protostar forming at the center, and showed the effects 
of an active accretion disk and an empty bipolar outflow 
cavity on the SEDs and the images. The stellar mass 
chosen is high enough that the star has a high luminos- 
ity (L = 2.8 X lO^L©) but low enough that the ionizing 
luminosity is small; we ignore the latter here. We in- 
cluded gas opacities inside the dust destruction front in 
the disk since hot gas becomes optically thick. How- 
ever, the outflow cavities in Paper I were left empty (an 
optically thin limit), which is not realistic since dust is 
expected to survive in these outflows; furthermore, un- 
der certain conditions of high mass outflow rates, the gas 
in the outflow may become optically thick near the star. 
In this paper, we shall study the effects on the SED of 
the gas and dust in the wind filling the outflow cavities. 
We shall also consider how the outflow removes angular 
momentum and accretion energy from the disk. In the 
next section, we discuss our model setup. We present 
the SED and image results of our models in Section [3l 
and summarize our main conclusions in Section IH In fu- 
ture papers, based on the fiducial model developed here, 
we shall present the SEDs and images of an evolutionary 
sequence of MYSOs, and also examine their dependence 
on the initial conditions of the core. 

2. MODELS 

We now briefly describe the model setup in our pre- 
vious work, which we continue to use here. For a more 
detailed description of the model we refer readers to Pa- 
per I. Following MT03, a star- forming core is deflned as 
a region of a molecular cloud that forms a single star or 
a close binary via gravitational collapse. The core is as- 
sumed to be spherical, self-gravitating and in near virial 
equilibrium. In the flducial case presented here, we study 
a core with an initial mass of = GOM©. The core 
is surrounded by a larger, self-gravitating clump with 



a mean surface density Hci = 1 g cm ^, which corre- 
sponds to pressure at the surface of the core of Pd = 
0.88GS2^ = 5.9 X 10-^(Sci/gcm-2)2 dyn cm-^. Such 

a core has a radius = 5.7 x 10~^(i;ci/g cm~^)~^/^ 
pc (1.2 X 10^ AU in the flducial case). We also con- 
sider two variants of this model in which the surface 
density is smaller or larger than the flducial value by 
a factor of 0.5 dex. We assume that the initial density 
distribution in the core is a power-law in spherical ra- 
dius, p (X r~^p ^ and adopt hp = 1.5. Such structures 
have recently been conflr med to be present i n Infrared 
Dark Clouds (IRDCs) bv lButler (|2QT2h . Here we 

study such a flducial core at the particular moment when 
the growing protostar reaches 8 Mq. Such a protostar 
has a radius = 12.05 i?© and a nuclear luminosity 
= 2.81 X 10^ Lo (MT03). 

The collapse of the core is described wit h an 

inside-out expansion wave solution (|Shul 119771 
iMcLaughlin fc Pudritzl I1997D . The expansion wave 
front moves outward, causing the material inside of it 
to collapse; the outer par t of th e core remains static. 
We use the solution by lUlrichI (|1976l hereafter the 
Ulrich solution) to describe the effects of rotation on 
the velocity fleld and the density proflle inside the 
sonic point, where the infall becomes supersonic. In 
this solution, angular momentum is conserved until the 
material falls onto the midplane and forms an accretion 
disk. Assuming the rotational energy of the core is 2% 
of its gravitational energy {/3 = ^rot/|^grav| =0.02), the 
radius of the disk Wd is 4.5 x 10^ AU in the flducial case 
(Eq. 13 in Paper I). 

Given a sufflciently large effective viscosity, the disk 
transfers angular momentum outwards and enables the 
high accretion rates required to form a massive star. The 
disk around a massive protostar can have a high mass. 
In all our models with disks, we assume the mass ratio 
between the disk and the protostar is a constant {fd = 
rrid/m^ = 1/3), motivated by an expected rise in effec- 
tive viscosi ty due to disk self-gravity at about t his value 
of fd^ (e.g. lAdams et al.][T989L IShu et al.l ll990l . IGammid 
I2OOII ). As in the previous paper, the di sk structure is de- 
scrib ed with an "a-disk" solution tShakura fc SunvaevI 
[T973h . but we have generalized this to include the effects 
of the outflow and the accretion infall to the disk (Section 

ED). 

The outflow sweeps up part of the envelope and thus 
controls the efflciency with which the star forms out of 
the core. We use the same basic setup as in Paper I, but 
we have reduced the total mass loading rate of the wind 
relative to the stellar accretion rate to = mw/m^ = 
0.1 (rather than 0.6 in Paper I) , which is a typ ical value 
for disk wind models ( Konigl fc Pudrit"3 l2000l also see 
Section [2^ . For better comparison to Paper I, we keep 
the opening angle of the outflow cavity to be 51°. These 
assumptions then lead to an increased star formation ef- 
flciency of e*/ = 0.64 (rather than 0.5 of Paper I). The 
major difference from Paper I, where we treated the out- 
flow cavity as being completely optically thin, is that we 
have developed a theoretical wind solution to evaluate a 
more realistic distribution of gas and dust in the outflow 
cavity, which will be discussed in Section [221 

2.1. Disk 
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In Paper I, we considered a standard Shakura-Sunyaev 
disk with a constant accretion rate and a zero torque 
inner boundary condition. However, we will now con- 
sider several additional effects that may be important in 
protostellar accretion disks: 1) the accretion rate varies 
with radius because material in the outflow (or the in- 
flow) leaves (or joins) the disk; 2) similarly, outflows and 
inflows alter the angular momentum distribution in the 
disk; 3) the disk itself is massive and contributes to the 
gravitational potential; and 4) the potential changes with 
time due to the growth of the disk and the central ob- 
ject. As in Paper I, we to neglect any transfer of angular 
momentum between the disk and the star via magnetic 
torques. Furthermore, we continue to ignore the effect of 
stellar irradiation on the structure of the disk, although 
we do include its effect on the emission by the disk. 

From conservation of the mass in an annulus at tu, we 
have 



d 

7 1 

dt dw 



dt 



(1) 



where ^(tu) is the surface density of the disk at w and v-^ 
is the radial velocity of the accretion flow (negative for 
inward motion). Here we have supplemented the stan- 
dard equation of mass conservation for disks (e.g., Frank 
et al. 2002) with w times dT^in{^)/dt, the mass loading 
rate per unit area of the disk from the (Ulrich) inflow, 
and w times dT^w{w) / dt^ the mass loss rate per unit area 
of the disk into the wind. Integration of this relation im- 
plies the accretion rate through the disk, rhacc, varies 
with radius w. 

macc(^) -2ttwY^v^ 

d 



dt 



md{nj) - min{ru) + mu,{ru)] + m*, (2) 



where dmd{vo) / dt is the rate of change of the mass of 
the disk inside tu, drriw/dt and dm[n/dt are the outflow 
and inflow rates inside tz7, and m* is the accretion rate 
onto the star. The disk extends from the stellar surface 
at lu* to We use the Ulrich solution to determine the 
rate at which mass flows from the envelope onto the disk. 
For consistency with Paper I, we adopt the same stellar 
accretion rate, m* = macc(^*) = 2.4 x 10~^Mq yr~^, for 
most of the models we consider. 

Similarly, we supplement the standard equation of an- 
gular momentum conservation for disks (Frank et al. 
2002) with terms representing the gain and loss of angu- 
lar momentum due to inflow and outflow, respectively. 



d 
dt 



d 



dvj 
1 dGT 

2iT dvo 



d^ 
' dt 



d^ 

' dt 



(3) 



where Gt = 2t:wvT,w'^ {dVt / dw) is the viscous torque, 
jin('^) is the specific angular momentum of the inflow 
joining the disk, jw{^) is the specific angular momen- 
tum (both rotational and magnetic) of the wind, and ft 
is the angular frequency. We have assumed that the disk 
extends to the surface of the protostar, and we further 
assume that there is no viscous transport of angular mo- 
mentum between the disk and protostar (i.e., Gt = 
at the inner edge of the disk). Combining Equations ^ 



and (|3|) we have 



dt 



dj_ 
dt 



1 



dj 
' dw 
dGT 



(jin -j) 9Sin , {j -jw) 



{27rwTi) dw 



dt 



dt 



-C4) 



where j = w'^Q is the specific angular momentum of the 
disk gas. We see that j decreases (angular momentum is 
transferred outwards) (1) due to the viscous torque, pro- 
vided dG^/dw < 0; (2) the material falls onto the disk 
with sub-Keplerian speed; and/or (3) the wind leaves 
the disk with a super-Keplerian speed due to magnetic 
torques. Note that the deepening of the gravitational 
potential well because of the growth of the star and the 
disk enters into dj /dt ex dVt/dt. 

As mentioned in the previous section, we assume the 
disk mass always a constant fraction of the protostellar 
mass. 



rrid rfidiwd) . 
Jd = — = \ = I/O, 



(5) 



where fadiwd) is the rate of change of the total mass in 
the disk. dT^w/dt and are given by the wind solution 
(see Section 12. 2p . with the total mass loss rate to the 
wind also assumed to be a constant fraction of the stellar 
accretion rate, 



(6) 



Then the total mass loading rate from the inflow to the 
disk is 

rn-.ni'^d) = (1 + + /d)^*, (7) 

with the mass loading rate per unit area (dT^in/dt) and 
its specific angular momentum jin given by the Ulrich 
solution. With these boundary conditions, from Equa- 
tions ([1]) and ([3]), ^(tu), dT^{w) / dt^ and other profiles of 
midplane temperature, disk scale height, can be solved 
with the Qf-disk model. A detailed formulation can be 
found in Appendix [A] Here we assume a is a constant 
throughout the disk that is self-consistently determined 
by these boundary conditions. In the final fiducial model 
of the series a = 1.37. 

Although we do not explicitly use an equation of global 
energy conservation in constructing our model, we note 
that it can be written as 



Gm*m* 



La 



■ Lw ^ Ld^ Edj 



(8) 



where G is Newton's gravitational constant, the l.h.s of 
this equation is the total potential energy released in 
the accretion process from infinity to the stellar surface. 
Lace = Gm*m*/(2r*) is the energy emitted as radiation 
when the accreting material hits the surface of the pro- 
tostar, is the total mechanical energy carried away 
by the wind, Ld is the energy radiated from the disk 
due to viscous dissipation (note that radiation from star 
or part of Lace may be absorbed and remitted from the 
disk; they are not counted in Ld)^ and Ed is the rate of 
change of the energy content (gravitational and kinetic) 
of the disk, which is actually much smaller (by 2 orders 
of magnitude) than L^j and Ld in our case. For the disk 
wind, Luj = l.22Ld in the fiducial model (see Section 
12. 2p , which means the disk luminosity is now smaller by 
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a factor of 0.45 compared with that in Paper I, since 
more than half the disk energy is extracted by the wind. 
The accretion luminosity Lace as it transitions from the 
disk to the star (2.45 x 10^ Lq in the fiducial case) is 
assumed to be emitted from the stellar surface homoge- 
neously, which indicates an effective surface temperature 
of the protostar of T^^acc = 1-43 x 10^ K, including the 
internal luminosity (2.82 x 10^ Lq) from the protostar. 

2.1.1. Thin Disk 

The disk profile is calculated numerically as described 
in Appendix [Al with the opacities used in our model. Un- 
like in Paper I, where we then fit the density and scale 
height with power-law functions of radius for implemen- 
tation in the radiative-transfer code, here we use the nu- 
merical solutions directly. The accretion luminosity is 
emitted from the disk, obeying the profile 47rT^/(3r) , 
where r is the optical depth from the midplane to the 
surface (see Appendix [A]). The midplane temperature, 
Tc, is calculated from the a-disk model. We find that 
it is also consistent with the temperature from the final 
radiative transfer simulation, since the heating from the 
star has negligible penetration to the disk midplane. 

The boundary between the infiow and the outfiow is 
described by an lUlrichI ()1976f ) streamline, which is the 
trajectory of an infalling particle with an initial angular 
momentum but with no pressure forces: 

^ cos ^-j^,^ esc sin Onj^QQQ . s 

COS ^^y, esc — COSO 

Here vud^^^^ ^w,esc is the radial distance to the point at 
which the trajectory crosses the disk plane {0 = 7r/2); 
the disk radius is determined by trajectories in the plane 
(^w,esc '7r/2), which gives r ^ Wd- Therefore, under 
the condition of a thin disk, the boundary between the 
wind launching region and the infiow joining region of 
the disk is zud sin^ O^^esc^ which is 270 AU in the fiducial 
case. 

Figure [T] compares the radial profiles of the thin disk 
in the various models. There are jumps in the profiles of 
density, scale height and temperature due to the drop of 
the opacity from dusty conditions at 1400 K to gas-only 
conditions at 1600 K and then a rapid increase of the 
gas opacity with temperature to ~ 10^ K. The material 
is piled up just before this transition area at T < 1600 K. 
The red dashed curves are numerically calculated with a 
standard a disk with a constant accretion rate. In P aper 
I (Model 8) and in Model 9 of this paper (see Sec. [23]), we 
used the fitted power laws, shown as the red dotted lines. 
The dashed blue curves (Model 10) show the effects of 
the radial variation in the accretion rate due to the inflow 
and the outflow, the effects of the angular momentum of 
the inflow (but not of the outflow), the effects of a time- 
dependent disk, and the effects of a non-Keplerian poten- 
tial. The flrst panel shows the accretion rate through the 
disk (Eq. [2]). Inside 100 AU, rhacc increases due to the 
mass loss to the wind {dm^/dt). From 100 AU to 270 AU 
the negative drrid/dt (see the discussion in Appendix [A]) 
becomes dominant, so rhacc starts to decrease. Outside 
270 AU, rhacc drops rapidly due to the inflow {dra{xi/dt). 
One effect of this spatially varying accretion rate is that 
the surface density of the inner disk (inside 200 AU) is 
less than that of the outer disk. The black curves show 



^ 2.8 




100 I 

0.1 1.0 10.0 100.0 

cj (AU) 



Fig. 1. — The accretion rate, density, scale height, surface den- 
sity, temperature profiles, and luminosity radiated from a thin disk 
inside cylindrical radius vu. riHe = O.lriH is assumed here. The 
dashed red curves are from a standard a disk with a constant ac- 
cretion rate. The dotted red curves are fitted power laws, which are 
used in Paper I (Model 8) and in Model 9 in this paper. The blue 
curves contain the effects of varying accretion rate in the disk due 
to the infiow and the outfiow, angular momentum of the infiow, 
disk growth, and a non-Keplerian potential (Model 10). Effects 
of angular momentum and kinetic luminosity of the outfiow are 
considered in Model 11 (black lines). 

the effects of including the wind torque. Panels 2-5 show 
that this brings the dust/gas transition inward, while the 
last panel shows that it reduces the viscous dissipation 
and the corresponding luminosity radiated from the disk 
by a factor of 0.45. 

2.1.2. Effect of Disk Thickness 
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The finite thickness of the disk has a substantial effect 
on where the boundary between the outflow and inflow 
intersects the disk surface, so we develop a simple ap- 
proximation to allow for this. The vertical distribution 
of the density in the disk (tu < zud) has a Gaussian pro- 
file determined by the temperature at the disk mid plane 
(see Eq. 16 in Paper I). In the absence of outflows, the 
height of the disk at vud was defined as the height of the 
innermost Ulrich streamline. This sets a volume den- 
sity at the disk surface. In Paper I, the disk surface at 
interior radii was set by having this same volume den- 
sity. In the presence of an outflow, the disk height in the 
wind-launching region is defined by the point at which 
the density in the disk equals the density at the base of 
the wind, which we have somewhat arbitrarily defined as 
the point where the wind velocity is 1% of the Keplerian 
velocity (see Appendix B2). For disks of finite thick- 
ness, the innermost Ulrich trajectory intersects the disk 
surface at tu^^max > '^d sin^ ^^^^esc- Our method of de- 
termining the disk height in the wind-launching region 
becomes inconsistent with the condition that the wind 
extend to the innermost Ulrich streamline if that stream- 
line would intersect our model disk beyond its outer edge 
at Now, we model the disk as being thin, so that 
the radiative flux in the disk is vertical, and correspond- 
ingly the outer boundary of the disk is a vertical wall 
at zud. The problematic case-which includes the fiducial 
case we consider-corresponds to having the outermost 
wind trajectory exit our model disk through the verti- 
cal wall instead of the disk surface. In fact, the disk is 
likely to be thinner near the outer edge, since radiation 
can escape radially as well as vertically. Correspond- 
ingly, the density at the base of the wind is likely to be 
larger than the value we estimate in Appendix B via a 
self-similar model. We adopt the simple approximation 
for the problematic case that the height of the disk at the 
outer radius, ^d,max, is the same as that of the outermost 
wind streamline at that radius. Since we do not adjust 
the wind density, this leads to a density discontinuity in 
the outer disk such that the disk density is somewhat 
higher than the density at the base of the wind. We 
assume that the density at the top of the disk remains 
constant as the radius decreases, until the density at the 
base of the wind becomes large enough to match it; inside 
that point, the density at the top of the disk is defined to 
be the density at the base of the wind. This approxima- 
tion has a relatively small effect: it needs to be applied 
only for w > 400 AU, i.e., a relatively small fraction of 
the disk, and it leads to a density discontinuity of only 
^ 20% between the disk and the wind in this region. 

2.2. Outflows 

Bipolar outflows are u biquitous around protostars (e.g. 
iKonigl <SrPudrital2000f ). These outflows (or winds: we 
will use these two terms interchangeably in this paper) 
are believed to be driven by the magnetic fields thread- 
ing the rotating disk (or the star) and powered by the 
accretion. Current wind models can be categorized into 
two classes: disk winds, in which the wi nd flows along the 
magn etic field lines threading the disk ( Konigl fc Pudritzl 
l2000f ), and X- winds, in which the wind is launched from 
a narrow interaction region between the stellar magnetic 
field and the inner edge of a truncated accretion disk 
(|Shu et al..„2000.) . Although the properties of these two 



winds in the near-disk region are quite different, their 
momentum distributions become si milar at regions far 
from the disk, withp^ ex {rsmOw)~'^ (i Matzner fc McKe^ 
Il999f ). In Paper I, this property, together with an as- 
sumed star formation efficiency from the core, was used 
to determine the opening angle of the outflow cavities, 
which however were left empty. Here, we now include a 
detailed density distribution for such outflows. 

In this paper we consider only disk winds. We base our 
analysis on the wind calculated by iBlandford fc Pavn^ 
(1982, we hereafter refer to t his a s the B P wind) and 
refined by Contopoulos fc Lovelacd (|1994 ), but we alter 
this so that the outermost streamline of the wind, at 
cylindrical radius tz7max(^), coincides with the innermost 
streamline of the accretion flow. In the Ulrich solution 
for the accretion flow, this streamline intersects the disk 
midplane at = tz7^ sin^ ^w,esc- We assume that there is 
a central cavity in the wind due possibly to t he poloidal 
field of the protostar, as in the X-wind model (|Shu et al.l 
1995), or to a line-driven wind from the protostar; the ra- 
dius of this cavity is Wc{z). For X- winds, one can show 
that za7c{z) is qualitatively similar to a BP streamline; 
since we are considering a BP-like wind here, we assume 
that vjc{z) does in fact have the shape of a BP stream- 
line. Wind streamlines begin at a cylindrical radius zjq 
and at a height Zdi'cuo). We assume that the innermost 
wind streamline emanates from the surface of the star, 
^co = BP winds have a density at the base of the 
outflow that varies as tu"^/^ and a velocity that scales 
as the Keplerian velocity, vk oc . As a result, the 

mass loss scales logarithmically with the radius of the 
outermost streamline, 

= 47r / zuopv^dzuo oc ln(c7max,o/^^'co), 

(10) 

where we have assumed that the disk is thin, so that 
the mass flux is determined by the z-component of the 
velocity. The density in the outflow at a height z' above 
the disk surface is given in Appendix [Bl in terms of the 
cylindrical radius of the footpoint, wq. 

As discussed in the last section, the outflow extracts 
angular momentum and accretion energy from the disk. 
The specific angu lar momentum o f a self-similar disk 
wind is (e.g. ^Konid fc Pudritzl l2000l) 

ju, Vtw\, (11) 

where wa is the cylindrical radius of the Alfven surface 
(where the wind becomes super- Alfvenic). The ratio be- 
tween the Alfven radius and the footpoint radius of a 
streamline, zua/^q-, can be consi dered as a co nstant and 
has typical value of ~ 3 (e.g. iKonigl fc Pudritzl |20T)Q[ ). 
For the BP wind we approximate here (see Appendix 
|B|) . wa/^q = 5.48. For such a wind, the specific energ y 
of the wind, 6^, is given by (IBlandford fc Pavnel[T982l) 

2.3. Dust and Gas Opacities 

Dust shapes the SED of the protostar by transforming 
visible and UV radiation into the infrared. In Paper I, 
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Fig. 2. — The density and temperature profiles for the fiducial model (Model 13) at different scales. riHe = O.lriH is assumed. The white 
contours divide the disk, the envelope and the outflow; all three regions are shown in the middle column. The black region in the right 
column is outside the core, and is assumed to be a vacuum. The density drops rapidly with z at the base of the wind as the wind accelerates 
and moves outward. The dotted lines show the streamlines of the disk wind. Each interval contains 10% of the wind mass loss. The short 
black contour in the lower left-hand corner of the left-most temperature plot is the dust destruction front (T = 1600 K). The dividing line 
between dusty and dust-free gas in that plot is apparent in the temperature jump at 10-15 AU, this boundary starts at about 2 AU. 



the dust was assumed to be completely destroyed by sub- 
limation at a temperature Tgub = 1600 K, which made 
the transition between dust and gas opacities relatively 
sharp. In order to improve the temperature convergence 
we make a smoother transition between the dust and gas 
opacities by linearly reducing the fraction of the dust 
content from 100% (gas-to-dust ratio 100) at 1400 K to 
0% (no dust content) at 1600 K, approximating the re- 
sults of Semenov et al. (2003). We identify 1600 K as a 
reference sublimation temperature. 

We assume that the wind is dusty if it is launched out- 
side the radius at which dust sublimates, neglecting other 
dust destruction processes such as shocks. Even though 
the temperature in the the wind launched from the in- 
ner dust-free disk soon drops below the dust sublimation 
temperature, this part of the wind is still likely to remain 
dust-free since dust grains are not expected to form in 
such a short time at these densities. We therefore define 
the boundary between the dust-free and dusty outfiow by 
a streamline starting from the disk surface where the sur- 



face temperature is at the sublimation temperature (1600 
K). In our fiducial model, this radius is ^2 AU. Figure 
[2] shows the streamlines (white dotted curves) and the 
density profile (color scale in the upper panels). Each in- 
terval of the streamline contains 10% of rhw The dusty 
and dust-free parts of the wind can be recognized in the 
temperature profile. About 66% of the wind is dusty. 
Note because of the finite disk thickness, the outermost 
wind streamline emerges from the disk surface at as 
discussed in Section 12.1.21 On the other hand, we used 
the thin-disk equations in calculating the dynamics and 
vertical profile of the disk, and in that approximation the 
wind emerges from the midplane at zud sin^ 0^ sec (see 

Fig. HD. 

We divide the dusty part of the simulation into four 
different regions: (1) the envelope, (2) the low density 
regions of the disk (nn < 2 x lO^^cm"^), (3) the high 
density regions of the disk, and (4) the dusty part of the 
outflow. Opacity models for the first three regions were 
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TABLE 1 

Properties of features included in the model series. 



= SMq, = 2.8 X 1O^L0, Mc = 60Mq, Rc = 0.06 pc, ^-a;,esc = 51°, vjg = 4-5 x 10^ AU 



Models 


Opacity 


Disk 


Outflow 


8 


gas yl > ioUU IV J 
dust (T < 1600 K) 


standard ol disk model 
with constant accretion rate 


empty 


9 


smoottiGr trcinsition 
(1400 K < T < 1600 K) 


as above 


empty 


10 


as above 


effects of inflow / massive growing disk on 
accretion rate / angular momentum distribution; 
effects of outflow on accretion rate 


empty 


11 


as above 


including effects of outflow in extracting angular 
momentum and accretion energy from disk 


empty 


12 


as above 


as above 


dust free 


13 


as above 


as above 


dust + gas 



TABLE 2 

Parameters of models comparing cores with different surface density. 





Model 131 


Model 13 


Model 13h 


initial core mass Mc {Mq) 


60 


60 


60 


mean surface density Ed (g/cm^) 


0.316 


1 


3.16 


core radius Rc (pc) 


0.10 


0.057 


0.032 


outer radius of disk vjd (AU) 


801.4 


449.4 


253.4 


stellar accretion rate m* {Mq /yr) 


1.035 X 10-4 


2.398 X 10-4 


5.667 X 10-4 


stellar radius r* {Rq) 


11.3 


12.0 


5.93 


isolated stellar surface temperature T* (K) 


1.25 X 10^ 


1.22 X 10^ 


1.74 X 104 


isolated stellar + accretion temperature T*,acc (K) 


1.37 X 10^ 


1.42 X 10^ 


2.62 X 104 


isolated stellar luminosity L* {Lq) 


2.82 X 10^ 


2.81 X 10^ 


2.84 X 10^ 


stellar accretion luminosity Lace ( Lq) 


1.13 X 10^ 


2.45 X 10^ 


1.18 X 104 


disk luminosity ( Lq) 


5.55 X 10^ 


1.10 X 10^ 


5.44 X 10^ 


total luminosity Lboi (^0) 


4.51 X 10^ 


6.46 X 10^ 


2.01 X 104 



described in Paper 1. The opacity model for the outflow 
has smaller grains than that in the envelope and has no 
ice mantles (Kim et al. 19 9J). These d ust opacity models 
are default options in Whitney et alls (2003a) radiation 
transfer code. 

Gas opacities are used in the dust-free, hot inner disk 
and in the outflow launched from it. The mean gas 
opacity at ~ 1500 K is about 3 orders of magnitude 
lower than the opacity of the dust at ~ 1400 K, but be- 
comes comparable or even higher when the temperature 
reaches ~ 10^ K, where the gas becomes collisionally ion- 
ized (note that we are not including photoionization from 
the protostellar radiation field) . For temperatures higher 
than ~ 3000 K, we adopt gas opacities from the OP 
Project (e.g. Seaton| [2Q05l). For regions with T < 3000 K, 
the opacity model is from I Alexa n der fc Ferguson ! ([r994). 
Howev er, the opacity given by [Alexander fc Ferguson 
{"1994") has a minimum at 2000 K and increases when 
T drops below that because of molecule and dust forma- 
tion. Since we expect minimal dust formation to occur 
in the rapidly expanding outflow, in its dust-free region 
where T < 2000 K we adopt the opacity at 2000 K. More 
details about these gas opacities are given in Paper I. 

2.4. Simulations 



to perform our calculations. To calculate the equilib- 
rium temperature, we use the algorithm bv iLucvl (| 19991 ), 
which is implemented in the new version of the code. It 
sums the pathlengths of all the photons passing through 
the cell rat her than counting only t hose that are ab- 
sorbed (e.g., iBiorkman fc WoodI [20011 used in Paper I), 
and it thus can quickly reach temperature convergence 
with fewer photon packets. Starting from a uniform and 
cold (T = 0.1 K) state, the temperature profile becomes 
stable within ~ 5 iterations with a number of photon 
packets equal to the number of grid cells (~ 1.5 x 10^). 
There are still some oscillations at the transition regions 
between gas-only and dusty opacities because of the U- 
shaped opacity curve as a function of temperature, but 
they occur mostly in small, localized regions close to the 
disk midplane. The temperature profile on the disk sur- 
face and in the outflow is quite stable. We iterate 20 
times, average the temperature profiles of the last 10 it- 
erations, and then use this averaged temperature profile 
as the equilibrium temperature for the final iteration. 

Corrections made by adiabatic (expansion) cooling and 
(compressional) heating in the outflow and in the accre- 
tion flow are also considered. We assume that the tem- 
perature of the gas varies only slowly, so that it is ap- 
proximately in thermal equilibrium. The thermal energy 



We use the latest versio n of the Monte Carlo radiation 
transfer code by .Whitnev et al.l (|2003bl 2012 in prep.) 
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equation for cehs in these regions is then 
AaT^np{T)p + PgasV • v + V • {uw) 

N V V 



(13) 



where /tp is the Planck mean opacity, v is the velocity 
field, Pgas is the pressure of the gas, u is the internal 
energy of the gas, L is the total luminosity, N is the to- 
tal number of photon packets, V is the volume of the 
cell, and p'^nyl sums all the optical depths that pho- 
ton packets have traveled across the cell (jLucvl 11999) . 
From left to right, the four terms in the equation are 
the emission, the adiabatic heating or cooling, advec- 
tion, and the absorption respectively. Adiabatic cooling 
becomes important in the outer region (r ~ R^) of the 
outflow, where the temperature drops to < 100 K, since 
the adiabatic cooling rate is proportional to T whereas 
the thermal emission term is proportional to T^. In these 
regions adiabatic expansion reduces the temperature by 
a factor of 2 to 3. Adiabatic heating in the infalling en- 
velope is unimportant (the emission term is ~ 10^ — 10^ 
larger than the adiabatic heating term). Note that Equa- 
tion (p!3]) is valid for both the gas and the dust only 
in regions where they are well coupled, which happens 
for densities nn > 10 ^ cm~^ in dense molecular clouds 
(|Dotv fc: NeufelJl997 ). This density threshold is at a ra- 
dius ~ 1000 AU in the outflow. Beyond this radius, the 
calculated value of the dust temperature should remain 
valid, although the gas will become colder. The tempera- 
ture profiles are shown in Figure O The low-temperature 
region in the outflow that extends to a cylindrical radius 
of - 100 AU at height of 600 AU is due to the lower 
opacity of the dust-free gas. The black region along the 
axis is the inner cavity due to the stellar magnetic field 
or the line-driven stellar wind. 

A 1000 X 1500 grid is used to resolve the azimuthally- 
symmetric rxO space. In r space, 750 cells are used to 
resolve the disk with a finer grid at smaller r (~ 400 cells 
are used to cover the region inside 5 AU). In space, an 
especially fine grid is made to resolve the midplane as 
well as the inner and the outer boundary of the out- 
flow. About 700 cells are used between 20° above and 
below the disk midplane. About 100 cells are used within 
l9 = 5°, and - 300 cells between 6 = 50° and 53° (the 
opening angle of the outflow cavity is 51°). For each 
run, a high S/N SED and images at different wavelengths 
and instrument bands {SpitzerlRAC and MIPS, 2MASS 
J, H, K bands, GTC CanariCam, Herschel PACS and 
SPIRE, and SOFIA FORCAST) are produced for view- 
ing angles of 30° and 60°. For most of the models shown 
here, 10^ photon packets are used to obtain SEDs and 
5 X 10^ photon packets are used to produce images. 

2.5. Model Series 

Following the same methodology of Paper I, we develop 
a series of models to show the effects of the improvements 
discussed above. We start this series with Model 8, the 
last model in Paper I, but now calculated with the new 
version of the code. In Model 9, we include the smoother 
transition between the gas opacities and dust opacities 
in the disk. In Model 10, we improve the disk by adding 
the effects of varying accretion rate in the disk due to the 
inflow and the outflow, angular momentum of the inflow. 
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Fig. 3. — SEDs for the model series, assuming a distance of 1 
kpc. Two typical inclinations of 60° and 30° are shown. The lower 
dotted blue curve is the input stellar spectrum (black-body), and 
the upper dotted one is the total luminosity of the star, including 
both the internal luminosity and the energy released by the gas as 
it accretes from the disk to the star. We assume that this radiation 
is emitted as a black-body. 



disk growth, and a non-Keplerian potential. The effect 
of the outflow in extracting angular momentum and ac- 
cretion energy from the disk is included in Model 11. In 
Model 12, we apply the wind solution to the outflow cav- 
ity but assume it to be dust-free (i.e., only gas opacities 
in the outflow cavity). We also start to consider the flow 
terms in the energy balance equation from here. Model 
13 is the fiducial model, which has all the elements dis- 
cussed above: a smooth transition between dusty and 
dust-free gas, an accretion disk structure that includes 
the effects of time dependence and outflows, and an out- 
flow that is dusty when launched outside the sublimation 
radius. The density and temperature structure of Model 
13 is shown in Figure (|2]). The properties of each model 
are summarized in Table [TJ 

We also construct two variants of Model 13, one with 
Sci = 0.316 g cm~^ (Model 131) and another one with 
Sci = 3.16 g cm-2 (Model 13h). The size of the core 
scales as S^"^^^, the density as S^/^' pressure as S^^, 

and the accretion rate as S^/^ (MT03). The parameters 
we use here correspond to a change in the clump surface 
density by a factor of 0.5 dex and the surface pressure 
by a factor of 1 dex. Table [2] compares the parameters 
of Model 131, 13, and 13h. 



3. RESULTS AND DISCUSSION 
3.1. SEDs 
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Fig. 4. — SEDs at an inclination of 60° and a distance of 1 kpc, 
comparing Model 12 and Model 13. The total flux and the flux 
from each component are shown. The two dotted blue curves are 
the input stellar spectra associated with the internal luminosity and 
the internal plus accretion luminosity, emitted as a black-body in 
each case. 

Figure [3] shows the SEDs of the models at inclination 
angles between the line of sight and the rotational axis 
of 60° and 30°. At an inclination of 60°, the line of sight 
passes through the dense envelope, which obscures the 
protostar and the disk. At an inclination of 30°, the 
line of sight passes through the outflow cavity; because 
the density in the cavity is lower than in the envelope, 
and because only part of the outflow is dusty, we can see 
much more deeply into the cavity than into the envelope. 
Some important features in the SEDs include the water 
ice feature at 3 //m and the silicate feature at 10 /im. As 
discussed in Paper I, the original output SEDs show some 
significant emission features such as the Ha and Paschen 
a lines, but we subtract these features and smooth the 
SEDs at wavelengths < 2/im, since the line profiles may 
not be exactly correct given the spectral resolution of the 
radiative transfer simulation and since this line emission 
is not the focus of our study. 

The differences of the SEDs among Models 8-12 are 
relatively small. The angle-dependent SEDs of Model 8 
are the same as in Paper I, showing that the new method 
of calculating the temperature produces the same results 
as the method used in the previous paper. Model 9 has a 
smoother transition between the gas and dust opacities, 
which is both physically more realistic and achieves bet- 
ter temperature convergence. Model 10 and Model 11 
update the disk profile by including the effects of the in- 
flow, the massive disk, and the outflow. Compared with 
Models 8 and 9, they have less dense but thicker inner 
disks, so that starlight can heat the inner parts of the 
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Fig. 5. — SEDs of the flducial model (Model 13) at four typical 
inclinations, assuming a distance of 1 kpc (F^ in the upper panel 
and uFiy in the lower panel). The thicker curves are total fluxes 
and the thinner red curves show only the scattered light. The two 
dotted blue curves are the input stellar spectra associated with 
the internal luminosity and the internal plus accretion luminosity, 
emitted as a black-body in each case. 



disk more easily, producing warmer disk surfaces, while 
the base of the envelope is more shielded and becomes 
cooler, thereby making the near- and mid-IR emission 
higher and the far-IR peaks lower. In Model 11, half 
of the accretion energy drives the outflow, reducing the 
emission compared to that in Model 10, particularly be- 
tween 1 fim and 10 /im. In Model 12, gas opacities are 
included in the outflow cavity. However, most of the out- 
flow is very cool so that the gas opacity is low, making 
the SEDs very similar to those of Model 11. 

The dust in the outflow cavity affects the SED more 
signiflcantly. As shown by the SED of Model 13 at 
a viewing angle of 30° in Figure O the emission at 
short wavelengths (A ^ 1/im) is now suppressed. The 
reprocessed emission comes out as extra near- and mid- 
infrared radiation. At larger inclinations, the fluxes at 
short wavelengths become higher than in the previous 
models due to scattering by dust grains. Figure 3] com- 
pares the SEDs of Models 12 and 13 at an inclination 
of 60° in more detail. Flux components of different ori- 
gins are also shown. The existence of dust in the outflow 
cavity makes the disk warmer, producing higher mid- 
and far-IR fluxes. The outflow also becomes hotter and 
emits much more infrared radiation, especially around 
10 to 20 /im, where the outflow emission becomes dom- 
inant. Depending on the sensitivity of the observations, 
some of the differences among these models may not be 
detectable, such as from Model 8-12. However, the dust 
in the disk wind (Model 13) can produce factors of sev- 
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eral differences to the SED in near- and mid-IR, which 
can be constrained by the observations. As we show be- 
low, dust in the outflow cavity also has significant effects 
on the infrared images, which are likely to be more use- 
ful for observationally testing these components of the 
theoretical model. 

The SEDs of the fiducial model (Model 13) at four 
inclinations are shown in Figure IsFl Here we extend 
the wind to 5i?c- If the wind extended to lOi^c, 
would increase by ~ 1% at an inclination of 10°, and 
~ 0.1% at 30°. Since the density of the wind drops as 
p oc tu"^ ~ (see Appendix [B]) . we conclude that the 
SEDs are not affected by truncating the wind at 5i?c- 
As a result of dust in the outflow, the optical emission 
is suppressed by a factor of 2-3 when the line of sight 
passes through the outflow cavity at a viewing angle of 
30°). However, at a viewing angle of 10°, the extinc- 
tion is negligible and the SED closely follows the stellar 
spectrum at short wavelengths. The SEDs (in vFj^) at 
these viewing angles are relatively flat from the near-IR 
to about 60 /im, which may be used as an identifler of 
near face-on massive protostars. The red curves show 
the SEDs in scattered light. If the line of sight passes 
through the envelope, all the light at wavelengths shorter 
than ~ 10/im has been scattered. Direct emission domi- 
nates for models with lines of sight passing through the 
outflow cavity, even at short wavelengths. 

As the extinction varies with inclination of the viewing 
angle, the IR (or bolometric) luminosity inferred from 
the observed flu x varies, which is known as the ^^flash- 
light effect" ( Na kano et al.lll995L lYorke Bodenheimerl 
[l999i) . Figure [6] shows how the inferred IR luminosities 
from the near-IR to the far-IR change with viewing an- 
gle. In the near-IR, the flux for a face-on view can be 
higher than that for an edge-on view by 3 orders of mag- 
nitude, and for wavelengths ^ 40 /im, it can be larger 
than the value averaged over the sky by a factor ~ 3. 
This contrast decreases as the wavelength increases; at 
100 /im and longer wavelengths, it becomes very small. 
When the line of sight passes through the outflow cavity 
close to face-on, the inferred IR luminosities (from the 
observed uFjy) in the (1-8) jam band give a good es- 
timate of the total bolometric luminosity, being low by 
a factor of ~2. Figure [3 shows how the bolometric lu- 
minosities inferred from the SED depend on the viewing 
angle and how they compare with the true bolometric 
luminosities. When looking down the outflow cavity, the 
bolometric luminosity inferred from the SED is higher 
than the true value by a factor of ~ 3 The inferred lu- 
minosity gradually decreases with increasing inclination 
angle and drops rapidly when the angle is large enough 

^ We note that in our highly ideahzed model, in which the axis 
of the outflow is exactly straight and the BP streamlines persist to 
large distances, it is possible to see the protostar through the dust- 
free part of the outflow only for viewing angles that are very small. 
For Model 13, in which the dusty outflow commences at vuq = 2 AU 
and the core size is Rc = 1.2 x IC^ AU, the central dust-free cavity 
has a size of about 600 AU at Rc, corresponding to an opening 
angle of about 3°. At greater distances from the protostar, this 
angle becomes even smaller. In fact, however, protostellar outflows 
are not observed to remain exactly straight to large distances, so it 
is quite possible that there is no inclination at which the protostar 
is completely unobscured by outflow dust. However, in the flducial 
model considered here, the amount of obscuration by outflow dust 
is negligible for moderate angles of inclination (~ 10°). 



that the line of sight passes through the envelope. At 
high inclinations, the inferred bolometric luminosity is 
lower than its true value by a factor of ~ 6. 

Figure [8] compares the SEDs of Model 13, Model 131 
and Model 13h, which have cores with the same mass 
but which are embedded in clumps of different surface 
densities, Sci- As noted above, the density in the core 

3/2 

scales as S^j , which has two effects on the SED: First, 
since the accretion rate scales inversely with the free- 
fall time, which varies as p~^l^ ^ the accretion luminosity, 

3 /4 

which scales with the accretion rate, varies as S^j ; this 
affects the SED at all wavelengths. Second, since the sur- 
face density of the core scales linearly with that of the 
clump in which it is embedded, the extinction due to the 
core is proportional to Sci; this affects the SED mainly 
at mid-IR and shorter wavelengths. Note, that we are 
not including the additional extinction that would result 
from the clump itself. So, as shown in the SEDs, Model 
13h generally has higher fluxes at all wavelengths than 
Model 13 and Model 131, due to the much higher luminos- 
ity. At lower inclination, the higher extinction due to the 
outflow makes the FUV fluxes in Models 13 lower than 
those in Model 131 (although these fluxes will in general 
be difficult to observe due to additional foreground ex- 
tinction from dust in the clump and the diffuse Galactic 
ISM). At higher inclination, the lower extinction of the 
core makes the mid-IR flux of Model 131 higher than that 
of Model 13. In this case, the mid-IR flux from the disk 
in Model 131 is higher and less buried by the flux from 
the envelope, producing a more prominent 10/im silicate 
absorption feature than in Model 13. This again implies 
that the depth of the silicate feature depends not only 
on the optical depth but also on the ratio of the mid-IR 
emission from the disk to that from the envelope. 

An additional factor is the dependence of the protostel- 
lar evolution (especially size of the protostar) on the ac- 
cretion rate. The protostar in Model 13h is much smaller 
(by a factor of about 0.5) in radius than the other mod- 
els: it has not yet reached the luminosity wave stage that 
swells the star (MT03). Thus, its accretion luminosity is 
signiflcantly higher than that expected from just its en- 
hanced accretion rate. The dependence of these results 
on the use a simplifled one zone protostellar evolution 
model (MT03) compared to a more detailed treatment, 
will be examined in a future paper. 

3.2. Images 

Figures [9lfT2l present the images for the flducial model 
at inclinations of 60° and 30°. We show both resolved 
images and those after convolving with the PSFs of 16 
different instrument bands. A distance of 1 kpc is as- 
sumed and no foreground extinction is applied. Here, 
5 X 10^ photons are used to produce images, but one still 
can see the noise due to the Monte Carlo method in the 
extended outflow, especially in the bands with low fluxes. 

These results conflrm the conclusion of our previous 
study: The outflow cavity has a signiflcant influence on 
the images up to 70/im, especially in the mid-IR, where 
the cavity wall (the inner part of the accretion flow that 
is heated by the protostar) dominates the morphology of 
the image. At around 20 /im, the thermal emission from 
the dusty outflow itself dominates the morphology. At an 
inclination of 60°, the receding side of the outflow is much 
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Fig. 6. — Inferred IR luminosities at different viewing angles and at different wavelengths (plus symbols and curves), compared with the 
true IR luminosities at these wavelengths (dashed lines) and the bolometric luminosities (dotted lines). The black plus symbols are at 20 
viewing angles with cos^view evenly distributed. The red plus symbols are at viewing angles of 0.1°, 1°, 5°, 10°, 30°, 60°, and 90°. The 
data shown with red plus symbols have a higher S/N than those in black. 



dimmer (by 1 to 2 orders of magnitude at 10 - 20 /im) 
than the facing side. The source becomes more symmet- 
ric at longer wavelengths. At ~ 40/im, the receding out- 
flow and cavity appear with a flux within a factor of 10 of 
that from the approaching side. The envelope dominates 
the image at wavelengths > 70/im. In addition to the 
wide angle outflow, a much narrower jet-like structure 
due to the denser region in the outflow close to the axis 
is also seen at an inclination of 60°. There is also a trend 
that from near-IR to ~ 20/im, the outflow appears to be 
narrower. Recall that we have assumed that the outflow 
extends to a radius of 5Rc. This extended outflow ma- 
terial is brightest in the near-IR because it is dominated 
by the light from the hot central region that is scattered 
towards the observer. Note that our model does not in- 
clude the material of the self-gravitating clump that is 
pressurizing the core. A real source would most likely be 
embedded in such a clump, and the appearance of the 



outflow outside the core would most likely be affected by 
additional foreground extinction. 

At an inclination of 30°, the star and disk could be 
seen directly in the previous models (up to Model 11) 
because the outflow cavity was treated as being optically 
thin. How are the locations of the r = 1 surfaces changed 
by the existence of the dust in the outflow cavity? Fig- 
ure [13] shows the r = 1 surfaces at different wavelengths 
seen from outside of the source at an inclination of 30° 
for Model 13. Even though there is dust in the outflow, 
the observer may still see down to the very base of the 
disk wind close to the protostar and disk even in the 
optical band. The r = 1 surface basically lies on the sur- 
face of the accretion disk in the near infrared bands, and 
reaches deeper as the wavelength increases. At inclina- 
tions large enough that the line of sight passes through 
the envelope, the source is relatively unobscured at wave- 
lengths ^ 100/im. With the angular resolution achieved 
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Fig. 7. — Inferred bolometric luminosity at different viewing an- 
gles (plus symbols and curves), compared with the true bolometric 
luminosity (the dotted line). The meaning of the black and red 
symbols are same as Fig. (|6]). 
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Fig. 8. SEDs of Model 131 (E^i = 0.316 g cm-^), Model 13 
(Eel = 1 g cm-2)^ and Model 13h (E^i = 3.16 g cm-^) at incli- 
nations of 60° (upper panel) and 30° (lower panel), assuming a 
distance of 1 kpc. The black curves are the total emission, and 
the blue curves are the input spectra of the protostar (internal + 
accretion). As discussed in the text, the protostar in Model 13h 
has a smaller radius, and therefore a larger accretion luminosity, 
than those in Models 13 and 131. 



by large sub-mm interferometers such as ALMA, such a 
disk (with a size of ~ 1000 AU) could be readily resolved 
at a distance of several kpc. 



Figure [14] shows the intensity profiles along the out- 
flow axis at inclinations of 60° and 30°. The profiles 
are convolved with a beam with FHWM of 4 " and nor- 
malized by the mean intensity at that wavelength of a 
2'' wide strip along the axis across the core. This strip 
extends to where the intensity drops to 10~^ of the max- 
imum intensity. As discussed above, at an inclination of 
60°, as one goes from shorter to longer wavelengths, the 
profile becomes more symmetric and the brightest point 
moves closer to the protostar. The extended facing out- 
flow {0 > 0) is brightest in the near-IR; the normalized 
brightness at 70 /im is similar to that at 10 /im. The op- 
posite side of the outflow can be best seen in the far-IR. 
At an inclination of 30°, the profiles are much more sym- 
metric. The degree of symmetry in the brightness profile 
appears to be a promising indicator of the inclination of 
observed massive protostars. 

Figure [T5l compares the images of Model 12 (dust-free 
outflow) and Model 13 (dusty outflow), showing that the 
dust makes both sides of the outflow cavities brighter at 
all the wavelengths we considered (which are between 2 
and 70 /im) except at 18 /im, where the emission of the 
dusty outflow dominates and the flux is more concen- 
trated to the central region. The dust in the disk wind 
also causes the outflow to appear narrower at 10 and 18 
/im than in the near-IR. Figure [16] shows similar inten- 
sity profiles as Figure [HI but perpendicular to the axis 
at an inclination of 60°. At 10 and 18 /im, the profiles 
of the bright component of the emission are narrower in 
Model 13. The difference between the two models be- 
comes small for A ^ 37/im. 

The effects of the surface density of the surrounding 
clump. Eel, are shown in Figure [17] The surface bright- 
nesses are normalized to the maximum value of Model 
13h at each wavelength. Model 13h has an envelope 
with higher extinction, which causes a higher contrast 
between the two sides of the outflow than in Model 13 or 
Model 131. The sides becomes more symmetric when 
is lower (Model 131). Also, we can see that the outflow 
cavity is brighter at mid-IR wavelengths (8 and 10.3 /im) 
in the two cases with higher extinction, as opposed to the 
model with the lowest extinction, for which the emission 
is concentrated at the center. 

4. SUMMARY 

We have constructed radiation transfer models for in- 
dividual massive star formation, with a 60 Mq initial 
core and an 8 Mq protostar forming at the center. 

1. Compared with Paper I, our model is improved in 
two important aspects: the effects of the outflow, inflow 
and time dependence on the disk surface-density pro- 
file are considered, and a detailed distribution of dust 
and gas in the outflow cavity has been studied in a self- 
consistent way. While most of the improvements in the 
disk model (except for the reduction in the accretion en- 
ergy, and hence the luminosity, due to driving the wind) 
have only minor influence on the results, the inclusion of 
dust in the outflow cavity can strongly affect both the 
SEDs and the images. 

2. We have presented a simple analytic model for the 
density distribution in a d isk wind that is based on the 
iBlandford fc Pavii3 (|1982l ) solution, generalized to allow 
for an arbitrary outer boundary. 

3. SEDs of the series of models at inclinations of 30° 
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10'^ 10'^ 10'" 10'^ 10'^ 10'' 10° 

Fig. 9. — Images for the fiducial model at an inclination of 60°, assuming a distance of 1 kpc, at different wavelengths. Each image has a 
field of view of 80^^ x 80^^ (80^^ corresponds to 8 x 10^ AU=0.4 pc), and is normalized to its maximum surface brightness, which is labeled 
in the bottom-left corner. The total fiuxes are labeled in the top-right corners. 
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Fig. 10. — Same as Figure [Q] except convolved with the beam of each instrument, labeled in the top-right and shown as red circles of 
radius equal to the HWHM on the bottom-right. The maximum surface brightness, which is dependent on the size of the beam, is labeled 
in the top-right corners. 
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Fig. 11. — Same as Figure (9] except at an inclination of 30°. 



and 60° are presented. Including dust in the outflow 
cavity significantly influences the SEDs at both incli- 
nations. The dust in the outflow suppresses emission 
at short wavelengths, while contributing a significant 
amount of near- and mid-IR emission. The disk also 
becomes warmer due to the re-radiated light from the 
outflow material, producing higher mid- and far-IR emis- 
sion. Close to a face-on view, at inclinations from ~ 10° 
to ~ 30°, the SED from near-IR to about 60 /im is rela- 



tively flat, which may be used to identify face-on MYSOs. 
The 10 /im silicate feature appears only if the mid-IR flux 
from the disk is high enough (e.g. due to lower extinc- 
tion of the envelope) so that it is not buried in the flux 
from the envelope, which indicates that the depth of this 
feature is determined both by the total optical depth at 
10/im and by the relative brightness of the disk and the 
envelope. 

4. We have presented images of the fiducial model 
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Fig. 12. — Same as Figure fTol except at an inclination of 30° 



for different observational bands, both resolved and con- 
volved with the resolution beams of the instruments. As 
discussed in Paper I, the outflow cavity is the most sig- 
nificant feature on the images up to 70 jam. The outflow 
cavity facing the observer may dominate the mid-IR flux 
and determine the morphology, while the opposite side 
can be seen only at wavelengths > 30 jam. Even though 
the outflow cavity is wide, at 10 to 20 /im, the dust dis- 
tribution in the disk wind can make the outflow appear 



narrower in the observed images. The distribution of 
the dust in the outflow affects the intensity profiles both 
along the axis and perpendicular to the axis. Comparison 
of such profiles with observation may potentially deter- 
mine the properties of the outflow. At a lower inclination 
(more face-on), the disk becomes visible even in the opti- 
cal bands, although foreground extinction from the natal 
clump and the larger-scale intervening ISM may limit de- 
tection of such emission in practice. 
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Fig. 13. — The solid contours show positions of the r = 1 surfaces for different wavelengths seen from outside. In this case the observer 
is looking from the upper right at an angle of 30° with respect to the axis. 





Fig. 14. — Intensity distribution along the projected outflow axis. 
The surface brightness are convolved with a beam of FMHW of 
4" and normalized by the mean value of a 2'' strip across the core 
along the axis at that wavelength. Different curves are for different 
wavelengths. The upper panel is at inclination of 60° and the lower 
panel is at 30° . The vertical dashed lines on the left and right mark 
the positions of an offset of half the core radius to the center. 



5. We studied the effects of the surface density (Sd and 
hence the ambient pressure, of the environment in which 
the core is embedded. Higher makes the core denser, 
leading to a higher accretion rate. It affects the SED 
in two ways: the higher accretion rate leads to higher 
luminosity, and the denser envelope (and outflow) leads 
to higher extinction, which can suppress emission at mid- 
IR and shorter wavelengths. The higher extinction also 
affects the images, making higher contrast between the 
two sides of the outflow, while in case of lower extinction, 
the two sides are more symmetric and the central region 
is more prominent. 

6. We studied how the inferred luminosity depends on 
the viewing angle (the flashlight effect) The luminosity 
inferred from the flux from a nearly face-on view can 
be higher than that from an edge-on view by 2 orders 
of magnitude in the near-IR and mid-IR. The inferred 
luminosity from a nearly face-on view is typically higher 
than the true IR luminosity by a factor of about 2. In 
the mid- and far-IR up to ~ 70 /im, the face-on fluxes, 
lyFjy^ can also give a good prediction of the bolometric 
luminosity. When integrated from a SED, the inferred 
bolometric luminosity is higher than the true value by a 
factor of 3 if looking down the outflow cavity, i.e. most 
of the radiation escapes from the polar direction. 
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Fig. 15. — Images at an inclination of 60° and a distance of 1 kpc, comparing Model 12 and Model 13. The size of the images is AO" x AO" . 
At each wavelength, the images are normalized to the maximum surface brightness of Model 13. The total flux and the maximum surface 
brightness of each image are also labeled in the upper-right and lower-left corners. 




Fig. 16. — Flux distribution along a strip crossing the core and perpendicular to the axis, comparing Model 12 and Model 13 at an 
inclination of 60°. Fluxes are convolved with a beam of FMHW of 4" and normalized to the values on the axis. From top to bottom, 
the intensity proflles are along strips crossing the core with offsets to the center of Rc/2 (the facing outflow), 0, and —Rc/2 (the receding 
outflow). From left to right, different wavelengths are shown. 
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APPENDIX 
A. AN IMPROVED a-DISK SOLUTION 

Here we summarize the formulae we use to include the effects of outflow/inflow and of a massive, grow ing disk in an 
Shakura-Sunyaev a-disk model. Following the approach in standard textbooks (e.g. iFrank et al.ll2002h . we write the 
conservation of the mass and angular momentum for an annulus at radius of vu as 



2.1 6 /zm 
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Fig. 17. Images of Model 131 (Ed = 0.316 g cm-^)^ Model 13 (Ed = 1 g cm-^)^ and Model 13h (Ed = 3.16 g cm-^) at an inclination 
of 60°, assuming a distance of 1 kpc. The size of the images is 60^^ x 60" . Each image is scaled to the maximum surface brightness of 
Model 13h at that band with their own maximum brightness labeled in the bottom left corners. The total fluxes are labeled in the top 
right corners. 



where H and dT^/dt are the surface density of the disk at w and its rate of growth, dT^^/dt and dT^in/dt are the 
surface mass loss/loading rates into the wind/onto the disk at tu, is the radial velocity of the accretion flow, 
Gt = 2t:wvT,w'^ {dVt / dw) is the viscous torque, and jin are the specific angular momenta of the wind and the 
inflow, and Vt is the angular frequency. Integration from gives 



^acc(^) 



d_ 

dt 



d 

m^{rln^) - — [Jd(tz7) + J^(tu) - JinM 



where subscript * indicates values on the surface of the star (i.e. the inner boundary of the disk), and 



dmd{w) 
dt 


/ ZTTtZJ ^— dU7 , 


dt 


[ 2ttw' ^-—^ dvo' ^ 


dJw{-^) 
dt ~ J 


1 2iTw jw dw , 


dJd{Tu) 
dt ~ J 


' 27iw' — {Y.w'^^)d^ 


We can write Equation 


(lA4l) as 



dt 



2i^w' — — ^ dw' ^ 
at 



dt 



I 27rw'ju 



dt 



■dw' ^ 



^acc(^) - rn^ 



37r- 



m* + rudiw) 



d [Jd{^) + Jw{^) - J\n{^)\ /dt 
^jG(rn^ + md{w))w 



(A3) 
(A4) 

(A5) 
(A6) 
(A7) 
(A8) 



(A9) 
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compared with a disk with constant accretion rate through R: 



The equation set for a Shakura-Sunyaev disk is: 



3^ 



1/2 



(AlO) 



(All) 



2.H = Cs/VL 

6.r = T^n{p,Tc) 
l.v = aCsH 

The radiation pressure term is smah for the current fiducial model, and thus neglected here. 

With the boundary conditions /^m*, md{^d) fdrn^, m^iwd) fw^*, rhini^d) = (1 + /d + fw)rn^, and 

dT^w/dt^ dT^in/dt^ jw and jin given by the solutions of the wind and infiow, at any moment t, on ce we hav e a profile 
of dT^ {w)/dt^ we can calculate T,{w) and other disk profiles (Tc, etc.) by combining Eq. ()A3p . Eq. ()A9p and Eqs. 
(|Alip . In fact, we set dT^/dt = at first, and calculate T^{w) at two close moments in time {t and t -\- dt)^ and use 
these two H profiles to re-estimate dT>{w)/dt, and iterate until the two converge. Note the size of the disk Wd can 
also vary at these two moments, i.e. the effect of an increasing disk size is also consistently included. We also use 
the surface density of the inner adjacent annulus to estimate S in the r.h.s of Equation ()A9p to make the calculation 
easier. 

The effect of the increasing size of the disk on the disk surface density can be estimated as below. The total derivative 
of the disk mass can be written as 



rnd{'^d)= / 27r — — — tz7(itz7 + 27rS(tz7^)tz7^d7^ 



dt 



dmd(w) 



dt 



■ 2TiTj(wd)wdWd' 



The size of the disk Wd oc md{vodY^^'^ (Eq. 13 in Paper I), so that Wd/vod = 1/(3 ■ 
we can write 



drridiru) 



dt 



^ rnd{wd) 



1 



3 kr, 



(A12) 

kp)[fad{wd)/md{wd)]. Then 
(A13) 



where S = m^i(tz7rf)/(7rtz7^). We can see that because the size of the disk is increasing, the mass of the disk inside 
Wd may either increase or decrease depending on the ratio of the surface density at the boundary of the disk to the 
average surface density of the disk. In our fiducial case, T,{wd)/^ ^ 1, then dmd{'^) / dt\^^^^ < 0, i.e. the disk mass 
inside Wd is decreasing. 

B. GENERAL FORMULATION OF THE DENSITY DISTRIBUTION OF THE OUTFLOW 

We wish to develop an approximate model for a disk wind that will fill the cavity described in Paper I. Protostellar 
disk winds are believed to be hydromagnetic in origin (Konigl & Pudritz 2000), but we shall not explicitly consider 
the magnetic field here. Our goal is to determine the dens ity distribution in the wind, and we shall make use of the 
results for the disk wind evaluated bv iBlandford fc Pa vne' (198 2l, BP wind) and refined bv iContopoulos fc Lovelac^ 
(l994). Our results asymptotically approach those of 'Ostriker (jl997l ) for self-similar cylindrical winds. We assume 
that the disk is thin and Keplerian, that the wind is steady and that it is launched from the surface of the disk in a 
self-similar fashion. 

We define the base of the wind as being at the point at which the vertical velocity Vz reaches some fixed, small 
fraction of the local Keplerian velocity vk] we label this velocity Vzq. The precise value of the ratio Vzo/vk is irrelevant, 
but the ratio must be sufficiently small that the base of the wind is close to the surface of the disk. We allow for the 
finite thickness of the disk, but we assume that it is small so that the velocity normal to the surface of the disk is 
approximately in the z direction. Let zu = r sin be the cylindrical radius of a streamline and wq denote the value of 
w at the base of the wind. The corresponding values of the innermost and outermost radii at the base of the wind 
are Wco and ti7max,0 7 respectively. The height of the base of the wind, i.e. the surface of the disk, is denoted by 
Zdi^o)- For a streamline starting from [tuo, Zd{vuo)]^ we define the effective height to be z'(mo ) = z — Zd(TUo). We 
assume that there is a central cavity in the wind created by the magnetic field threading the star (iShu et al lll995f ): let 
Wc[z\vUco)] be the radius of this cavity. We further assume that the wind is confined by the ambient medium, which 
sets tz7max[^'(^max,o)], the radius of the outer boundary of the wind; we shall specifiy this below. We introduce the 
normalized radius at the base of the wind. 
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We then have vk = vkc^q^^'^ and Vzo = VzcqXq^^'^^ where vkc = {Gm^/vucoy^'^ and VzcO are measured on the 
innermost streamhne. Let po('^o) be the density at the base of the wind for the streamhne originating at vuq. We 
assume that disk wind is launched in a self-similar manner, so that po(^o) varies as a power-law in radius, 

poi^o) = pco ( —] = PcoXq"^ , (B2) 

where pco is the density at the base of the wind along the innermost streamline. 

The mass-loss rate from a ring of width dzuco on both sides of the surface of the disk is 

drn^ =4:7TZAjopoVzodwo, (B3) 

= 47TwlQPcOVzcOX§~'^dXo , (B4) 

where we have made use of assumed power-law scalings of the density and velocity. The total mass-loss rate from the 
disk is then 

(B5) 



where 

'^max, 



In X] 



max, 



(B6) 



BP winds have ^ = f , and in this case /^^ = 1. 

The mass loss rate through rings of width dw at a height above and below the surface of the disk is approximately 
the same as that in the plane of the disk, 

drhw = AiiwpVzdw. (B7) 



Let 



n-o,.')-^^^^ (B8) 



be the ratio of the vertical velocity to the Keplerian velocity at the footpoint, and let 

9 In tz7o 



dhiw 



(B9) 



measure the change in the width of the ring. Equating the mass loss rate at a height z' to that in Equation (|B4[) and 
using Equation ()B5[) to eliminate pcO, we then find that the density in the wind is 

^max, 

This is a generalization of the results of lOstrikerl (|1997[ ) for cylindrical winds, which have rj = 1. Note that if the 
variation of xq is small (compared to that of txj), then Equation (IBlOp implies that p (x w' -2 (|Matzner fc McKedfTOOOl ) . 
In particular, this is the case for X- winds ( IShu et al.lfToM [19951 ). for which xqc^I. 

Bl. Approximate Streamlines 

In order to complete our estimate of the density in the wind, we need to specify the streamlines. We describe the 
streamlines in terms of the dimensionless radius, 

i?(t.o,.')-^^^, (Bll) 

and the dimensionless height, 

z' 

Z{wQ,z') = — . (B12) 

The innermost and outermost streamlines are Rc{Z) = Wc{Z)/wco and i?max(^) = '^max(^)/'^max,07 respectively. We 
assume that the streamlines emanating from points in the disk are a weighted mean of the innermost and outermost 
streamlines: 

tu(tuo, Z) = Wc{ZY-^w^UZY' (B13) 
Since this relation applies at the surface of the disk (Z = 0), one finds that 
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Thus, for the innermost streamhne {xq = 1), Equation (|B13p gives w = tUc, and for the outermost one {xq = Xmax,o) 
it gives m = tu^ax, as it should. The ansatzwe have adopted thus provides a smooth transition from the innermost to 
the outermost streamhne as the footpoint moves from Wco to tz7max,o- In terms of the dimensionless radius, we have 

R{Z) = R,{Zy-'R^,^{ZY. (B15) 

It is important to note that in both this expression and in Equation ()B13p . the radius of the innermost and outermost 
streamhnes are evaluated at the same value of the normalized height, Z, not at the same height above the disk, z\ 
In general, define the characteristic wind density as 

pch = -. — 2 — -ri ' (^^^) 

AtTW^qVkJw In Xmax, 



4.3 X 10-" /lO^y/YMoV^V ah, ] gem-3 (B17) 



Along a given streamline, the wind density (Eq. IB10|I is then 

p{wo,z') = pch 



ViZ)RHxo,z')xl 



(B18) 



The quantity 77, which measures the rate at which the streamlines separate, is given by Equation (jB9]), which can be 
expressed as 

\ dS 

^ ^ dlnZ dlnZ lnxmax,0 ' 

where all the partial derivatives are evaluated at constant z' . Note that if the streamlines are self-similar, so that R 
depends only on Z, and as a result Rc{Z) = i^max(^), this reduces to 

r]~^ = 1 — ^^1^^ (self-similar winds) (B21) 
am Z 

(Contopoulos & Lovelace 1994). 

B2. BP-like Winds 

We now focus on winds that are similar to those considered by BP. Such winds emanate from a disk with a density 
distribution that varies as with ^ = |, and as a result I^, = I (Eq. IB6p . For the particular numerical example 
BP considered, the dimensionless radius is given b}13 

Rc^l^U ln(l + 0.07Z), (B22) 

which agrees with the numerical values given bv iContopoulos fc Lovelac^ (|1994f ) to within about 10% for Z < 1246, 
beyond which the wind begins to contract toward the axis. For larger values of Z, this approximation indicates that 
the fiow is highly collimated, which is c onsistent with observations of jets from low- mass protostars (McKee fc Ostrikerl 
HoOTi). The dimensionless velocity ()B8|) in this wind is given by 

V{Z) :^ln(1.01 + 5Z0-^), (B23) 

which agrees with "Contopoulos fc Lovelace"' s (1994) numerical values to within about 10%, except for very small Z; 
we have inserted 1.01 instead of 1 in this expression in order to avoid the singularity that would otherwi se oc cur in 
the density at Z = 0. If both the innermost and outermost streamlines are BP streamlines, then Equation (jBlSp gives 
an approximation for the density in a BP wind. 

In reality it is likely that the outermost streamline is determined by boundary conditions, and we can use the analysis 
above to estimate the density in the resulting non-self-similar wind. As in Paper I, we a ssume that the outermost 
streamline of the wind corresponds to the innermost streamline of the accretion flow. ' UlrichI fl976) developed a simple 
model for a rotating infall characterized by a disk radius Wd- Gas originating far from the protostar at a polar angle 
^0 follows a trajectory in polar coordinates (r, 0) given by 

j_ ^ Ml -1^1)^ (J324) 

^ An expression that is more accurate at small Z is ~ 1 + rr^i u • ^ j-^- • • ^ 1 u- • n 1 + 

4 671nri+0 36Z+3 4x10-4^3)- in particular this expression gives physical conditions m the wind-launching region are likely to 

' \l ' ^ r / inr ^ 7^ i l^ • ^• ^ ^ ^ bc much morc compllcatcd than givcn by thc BP solutlou, howcvcr, 

the correct value or aR/aZ at Z = 0, where the wind is launched. j-u- ^^-j.- i • i ^ 

' ' so this additional accuracy is only apparent. 
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where /i = cos 6 and /jjQ — cos^o- Let Oyj esc be the largest angle at which the wind can escape, evaluated far from 
the disk. If this outermost streamline of the wind starts at (ti7max,0 5 ^d,max) and corresponds to a streamline of the 
accretion flow, then = ^w,esc^ and Zd^ma^ can be determined by 



"^max, 

One can then show that 



- ^of + ^d,max(l - Mo)/Mo + 2WdZd,maA^ - /io)/MO 
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(B25) 



(B26) 



where Z^ax = ^7^max,o = {z - 2;rf,max)/tx7max,o. Generally, ti7max,o can be values between Wd{l - /ip) a nd W d (see 
discussions in Sec. 12. ip . but in the fiducial model we present here, ti7max,o='^(i5 in which case, Eq. (|B25p and (jB26p 
become 
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1 



(B27) 



(B28) 



For O^^Qsc — 60°, we find that 77 is typically between 1 and 3, but can reach values as high as 7 close to the out- 
ermost streamline, for the case co nsidered in the text in which Xmax,o is very large. Fully expanded disk winds are 
asymptotically similar to X- winds (iShu et al.lll994[ ) and have 77 <C 1; such winds are possible only if Xmax,o is not large. 
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